function nominalbondprice_components = PriceNomBonds_Model20220926(parms, gesol, inflation_process, Tmax)
% see notes, inflation process is
% inf' = mu_pi + rho_pi*inf + F(N,lambda,z,z') + theta_pi*eps + eps', eps'~N(0,sigma_pi^2)
% P_nominal(N,lambda,z,X,Y,eps) = P_A(N,lambda,z)*P_B(X,eps)*P_C(Y)

    NN = parms.NN;
    Nx = parms.Nx;
    Nz = parms.Nz;
    
    % grid for interp2
    [NN_mesh,xx_mesh] = meshgrid(parms.grid_N(:),parms.grid_x(:));

    mu_pi = inflation_process.mu_pi;
    rho_pi = inflation_process.rho_pi;
    sigma_pi = inflation_process.sigma_pi;
    theta_pi = inflation_process.theta_pi;
    F_Nxzz = inflation_process.F_Nxzz; % F(N,x,z,z')
    
    % law of motion
    x_next_Nxzz = gen_law_of_motion_x(parms);
    Nnext_Nxz = gen_law_of_motion_N(parms,gesol);

    % Arrow-Debreau price(N,lambda,z,z')
    ADprices = permute(repmat(parms.StateTransitionProbs,[1 1 NN Nx]),[3 4 1 2]).*gesol.spd1period;


    P_A = zeros(NN, Nx, Nz, Tmax+1);
    P_A(:,:,:,1) = 1; % initialize

    P_A_next = ones(NN, Nx, Nz, Nz);

    inflation_factor = zeros(NN, Nx, Nz, Nz);
    
    as = zeros(1,Tmax+1);
    bs = zeros(1,Tmax+1);
    cs = zeros(1,Tmax+1);
    
    % initialize at mat=0
    as(1) = 0;
    bs(1) = 0;
    cs(1) = 0;

    for TTM=2:Tmax+1   % time to maturity    (n in notes)

        % P_A
        % Fprefactor = rho_pi*(1-rho_pi^(TTM-1))/(1-rho_pi); % INCORRECT? FORGOT SOMETHING???
        Fprefactor = 1 + rho_pi*(1-rho_pi^(TTM-1))/(1-rho_pi);
        inflation_factor(:) = exp(-Fprefactor*F_Nxzz); 

        for iznext = 1:Nz
            P_A_next(:,:,:,iznext) = reshape(interp2(NN_mesh,xx_mesh,P_A(:,:,iznext,TTM-1)',...
                Nnext_Nxz(:),x_next_Nxzz(:,iznext)),[NN Nx Nz]);
        end
        P_A(:,:,:,TTM) = sum(ADprices.*inflation_factor.*P_A_next,4);

        % P_B = exp(-a - b*X)
        as(TTM) = as(TTM-1) + mu_pi*(1+bs(TTM-1))- (1/2)*((1+bs(TTM-1)+cs(TTM-1))*sigma_pi)^2;
        bs(TTM) = (1+bs(TTM-1))*rho_pi;
        cs(TTM) = (1+bs(TTM-1))*theta_pi;

    end
    
    % save

    nominalbondprice_components.maturities = 1:Tmax;
    nominalbondprice_components.P_A = P_A(:,:,:,2:end);
    nominalbondprice_components.P_B_a = as(2:end);
    nominalbondprice_components.P_B_b = bs(2:end);
    nominalbondprice_components.P_B_c = cs(2:end);
    nominalbondprice_components.P_C_a = rho_pi*(1 - rho_pi.^(1:Tmax))/(1-rho_pi);
    nominalbondprice_components.P_A_format = 'P_A=P_A(N,x,z,T), T=1,2,...';
    nominalbondprice_components.P_B_format = 'P_B(X,eps,T)=exp(-P_B_a(T)-P_B_b(T)*X-P_B_c(T)*eps)';
    nominalbondprice_components.P_C_format = 'P_C(Y,T)=exp(-P_C_a(T)*Y)';

end